Decrease in air-sea CO2 fluxes caused by persistent marine heatwaves

Regional processes play a key role in the global carbon budget. Major ocean CO2 uptake at mid-latitudes counteracts CO2 release in the tropics, which is modulated by episodes of marine heatwaves. Yet, we lack essential knowledge on persistent marine heatwaves, and their effect on the CO2 sensitive areas. Here we show, using a 1985–2017 joint analysis of reconstructions, ocean reanalysis and in situ and satellite data, that persistent marine heatwaves occur in major CO2 uptake and release areas. Average air-sea CO2 flux density changes from persistent marine heatwaves are strongest in the Pacific Ocean with a 40 ± 9% reduction in CO2 release in the tropics linked to ENSO, and a reduction in CO2 uptake of 29 ± 11% in the North Pacific over the study period. These results provide new insights into the interplay of extreme variability and a critical regulating ocean ecosystem service, and pave the way for future investigations on its evolution under climate change.

Extreme events associated with unusually high-water temperature are ubiquitous in the global ocean. They can last from weeks to years, span from local to interbasin scale, and can reach depths of several hundreds of meters [1][2][3] . These so-called marine heatwaves (MHWs) occur due to either coupled air-sea interactions [4][5][6][7] , ocean internal processes such as horizontal and/or vertical circulation changes 8 , and are sometimes linked to large climate modes such as the El Niño-Southern Oscillation (ENSO) 5 . Over the past 35 years, MHWs have become longer-lasting, more intense and more extensive 9-11 very likely due to long-term anthropogenic change 1,9,10,12,13 .
Intense and long-lasting MHWs have been reported at different locations in the global ocean 3,5 . These include the 2013/2015 Northeast Pacific 'warm blob' 4,14 , the 1997/1998 El Niño 15 , the 2015/ 2016 Tasman Sea 16 or the 2012 Northwest Atlantic 17,18 . The duration of these major events ranges from several months up to 2 years, they are associated with dramatic increase in sea surface temperature, that can exceed 5°C in anomalies 19 , and they can extend over large regions, reaching sometimes~10 M km 2 5 . Due to their extreme nature, MHWs, and in particular the intense, persistent ones, pose a fundamental challenge for societies as they have devastating impacts on marine ecosystem and their services 1,11,20 .
The ocean acts as a net sink for atmospheric CO 2 and is absorbing about a ¼ of CO 2 anthropogenic emissions 21 (2.6 ± 0.6 PgC/yr over the 2009-2018 period), thereby mitigating global warming. CO 2 entering the ocean is then redistributed horizontally over large distances and into deep ocean layers where it is then stored for long time scales [22][23][24] . The magnitude and direction of the air-sea CO 2 flux density (F CO2 ) vary widely in space and time, and depend on hydrographic conditions, the ocean circulation system, biological net production and air-sea interactions. As a result, major CO 2 uptake areas are located at mid-latitudes, whereas CO 2 release takes place predominantly in upwelling areas such as the tropical ocean [25][26][27] (Fig. 1a).
Persistent MHWs linked to ENSO 5 affect the Tropical Pacific CO 2 source region and lead to significant reduction in CO 2 outgassing [28][29][30][31][32] . However, we lack essential knowledge about how major MHWs events affect other oceanic CO 2 sources and sink regions. Here we investigate the interplay and the impact of intense and long-lasting MHWs on the air-sea CO 2 flux density at the global scale. The study is built on a combined use of reconstructions from 1985 to 2017, direct measurements, remote sensing data and an ocean reanalysis. We first present the regions where particularly intense and long-lasting MHWs most frequently occur. We then quantify the impact of these extreme ocean events on oceanic CO 2 sink and sources areas. We further examine the interaction between these extreme ocean events and one critical oceanic CO 2 sink region in the North Pacific Ocean. Finally, we discuss these results with existing knowledge on the mechanisms in the Tropical Pacific region to obtain a large-scale view of the prevailing mechanisms driving coupled changes between one of the regulating ocean ecosystem service and extreme variability.

Results
Persistent marine heatwaves occurrence and oceanic CO 2 source and sink areas In the Tropical Pacific, intense and long-lasting MHWs (hereinafter denoted persistent marine heatwaves, PMHWs) have a strong impact on air-sea CO 2 fluxes 31 . We propose a specific new criteria to identify such PMHWs at the global scale based on the duration and the mean Sea Surface Temperature (SST) anomaly during a MHW. We first detect all MHWs that occurred from 1985 to 2017 by applying a standard MHW detection algorithm 2 to NOAA gridded SST data derived from AVHRR sensor 33,34 (see method section). The detection algorithm provides several metrics that describe MHWs, including the duration and the mean SST anomaly. Using these two metrics, we define PMHWs as MHWs whose duration and mean SST anomaly are greater than the 95th percentile of their global historical distribution, i.e., duration > 38 days and mean SST anomaly > 2.3 degrees Celsius. Finally, we focus on the regions where PMHWs have appeared several times over the past three decades, and that represent a recurring event effecting the ocean CO 2 sink, similar to El Niño events in the Tropical Pacific. To do so, we only consider the points where PMHWs have re-occurred at least three times during the 1985-2017 period (gray points in Fig. 1a)-which correspond to 25% of all grid points that have experienced at least one PMHW.
PMHWs most frequently occur in the largest oceanic CO 2 source and sink areas in the near-global ocean. The analysis is restricted to within 60°S and 60°N to exclude polar regions where periods of ice longer than 5 days obstruct the MHW detection algorithm (see method section). Critical sink regions are located at mid-latitudes in the Northern and Southern Hemispheres (solid contours in Fig. 1a), while ocean CO 2 outgassing predominantly occurs in upwelling regions such as the Tropical Pacific (dashed contours in Fig. 1a) and correspond to regions where climatological F CO2 is lower/greater (sinks/sources) than −1/1 molC/m 2 /year (solid and dashed contours in Fig. 1a Fig. 1a) are mainly located in the strongest oceanic sources and sinks areas in the near-global ocean, particularly in the Pacific Ocean ( Fig. 1a and Table 1). In the North Pacific and in the Tropical Pacific, the area impacted by PMHWs covers 3.8 and 2.9 * 10 6 km 2 , respectively, which corresponds to 21% and 16% of the total area of these basins. Whereas, in the North Atlantic and in the mid-high latitude Southern Oceans, the  and its Taylor decomposition (vertical bars). The contribution of temperature, Dissolved Inorganic Carbon (DIC), Alkalinity (ALK), salinity, wind and atmospheric partial pressure of CO 2 to F CO2 anomalies observed during PMHWs in the North Pacific CO 2 sink region for the 2009-2017 period were calculated using a first order Taylor expansion derived from the biogeochemical reanalysis (see methods section). The "total" bar corresponds to the sum of all contributing terms and corresponds to the Taylor approximation of F CO2 anomalies (black dot). The good agreement between the two implies that F CO2 anomalies are well approximated by the Taylor decomposition. The error bars correspond to the 95% confidence interval. b Contribution of horizontal, vertical diffusion, vertical advection, dilution and concentration due to freshwater fluxes, air-sea CO 2 flux density, a residual term and biological activity to the rate of change (tendency or trend) of DIC anomalies during PMHWs in the North Pacific CO 2 sink region for the 2009-2017 period (see methods section). The vertical bars represent the slope from linearly regressing each forcing term to the DIC anomalies trend 28 . A linear regression slope close to 1 indicates that a particular term produces in-phase anomalies of comparable magnitude. A slope near zero indicates that the term is not important in generating anomalies. The error bars correspond to 95% confidence intervals. area impacted by PMHWs covers only 10% and 5% of the total area of these basins.
The PMHWs impact on F CO2 is quantified using an ensemble of four observation-based products of F CO2 from 1985 to 2017 (see method section). We use a common approach using the ensemble spread as first order uncertainty estimate 37 . We would like to note though, that regionally, the uncertainty might be larger resulting from the lack of direct measurements. Previous studies, however, show that the uncertainty as a result of data sparsity is in the order of the ensemble spread adopted here 38 . Additionally data products are limited by 1×1 degree (2.5×2 degree for the JENA product) resolution by design and we are unable to formally quantify the role of the resolution on the uncertainty of our study. Amongst the largest oceanic CO 2 sinks and sources where PMHWs most frequently occur, the North Pacific and the Tropical Pacific are the most impacted. During PMHWs events, both areas suffer from a significant reduction in the air-sea CO 2 flux density with a change in uptake (29 ± 11%) and outgassing (40 ± 9%) respectively. In contrast, in the North Atlantic and the mid-high latitude Southern Oceans CO 2 sinks, the impact of PMHWs on F CO2 is close to 0 and negligible over the study period. We verified that, in these two regions, the small value is a general feature and is not just a cancellation of large positive and negative anomalies. The impact of PMHWs on F CO2 appears thus to be the most important in the Tropical and North Pacific, which is of considerable concern given their contributions to the global ocean carbon cycle [25][26][27] . The Tropical (dashed contour) and North Pacific (solid contour) represent a source and a sink of 0.52 and −0.60 PgC yr −1 , respectively which corresponds to −36 and 42% of the annual near-global net flux while covering roughly 5% of the near-global ocean area.

Persistent marine heatwaves and the North Pacific CO 2 sink
We use a state-of-the-art ocean biogeochemical reanalysis 39  the SST grid points that have experienced at least 3 PMHWS from 1985 to 2017. The exchange of CO 2 between the ocean and the atmosphere is driven by six variables 42 : wind, upper-ocean temperature, salinity, dissolved inorganic carbon (DIC) and alkalinity (ALK) as well as the atmospheric partial pressure of CO 2 . Temperature, salinity, DIC and ALK are influenced by hydrographic conditions, the ocean circulation system, and air-sea interactions. In addition, DIC and ALK are also influenced by the biological net production 42 . As observation-based products of F CO2 do not include these variables, the reanalysis becomes essential to pursue the analysis. The BGC reanalysis combines ocean circulation and biogeochemistry models together with in situ and satellite observations to provide a high degree of bio-physical realism 43 . The reanalysis skill is validated against the ensemble of observation-based products over the 2009-2017 period and in situ observations from an array of BGC-Argo floats during the 2013/2015 'warm blob' PMHW (see supplementary Information). The reanalysis shows good agreement with the observation-based products in estimating F CO2 anomalies due to PMHWs in the North Pacific. The reanalysis also agrees well with the float observations in reproducing anomalies in the four oceanic drivers known to control F CO2 (temperature, salinity, DIC and ALK) during the 'warm blob'.
During PMHW events, the reduction in F CO2 is the result of higherthan-usual temperature and negative DIC anomalies. We calculate a first-order Taylor series expansion of F CO2 anomalies to determine the contribution of the four oceanic drivers 28,44,45 (see method section). The Taylor decomposition (Fig. 2a) reveals that the reduced uptake of CO 2 during PMHWs in the North Pacific mainly result from the contribution of temperature (1.43 ± 0.02 molC/m 2 /yr), DIC anomalies (−0.81 ± 0.01 molC/m 2 /yr) and to a lesser extent ALK anomalies (−0.23 ± 0.01 molC/ m 2 /yr). The contribution from salinity, wind and the atmospheric partial pressure of CO 2 anomalies are small, and can be considered negligible. Sea surface warming during PMHWs reduces the solubility of CO 2 in the ocean resulting in a reduced uptake of CO 2 . In contrast, the decrease in DIC associated with a small increase in ALK (Fig. S3) enhances the uptake of CO 2 and as such offsets the thermal effect to the extent that the final F CO2 anomaly is~4 times smaller than it would have been from the thermal effect alone.
Next, we investigate the mechanisms leading to negative DIC anomalies. We examine the processes that drive the rate of change (tendency or trend) of DIC anomalies during PMHWs over the 2009-2017 period. The budget (or forcing) terms in the DIC trend equation consist of: horizontal and vertical advection, vertical diffusion, the air-sea CO 2 flux density, biological activity, dilution and concentration due to freshwater fluxes and a residual term (see method section). To highlight the contribution from each process to the DIC anomalies trend, we follow the method of Doney et al. 28 and examine the slope through linear regression of each forcing term to the DIC anomalies trends (Fig. 2b) (we verify that the intercept is approximately 0 because the average of the forcing term anomalies is null). A slope close to 1 indicates that a particular forcing term produces in-phase anomalies of comparable magnitude to DIC anomalies trend. In contrast, a slope near zero indicates that the term is not important, and a negative slope that the term produces out of phase anomalies.
Horizontal advection is the main driver for DIC anomalies. The linear regression slope of horizontal advection on DIC anomalies trends is the largest (0.70 ± 0.02 (unitless)) whereas the slopes of the other forcing terms are much smaller (<0.16 (unitless) for vertical diffusion anomalies and lower than 0.06 (unitless) for all the other terms). Furthermore and consistently with Ayers and Lozier 46 and Gruber et al. 23 , the reanalysis shows that, on average, there is a net horizontal divergence of DIC in the North Pacific CO 2 sink region (data not shown). The reanalysis suggests that the lateral removal of DIC is further accentuated during PMHWs, causing a decrease in DIC. In the North Pacific, extreme MHWs are associated with changes in horizontal advection due to wind speed modification 4,19 . This suggests that similar processes could potentially both drive thermal and DIC changes during PMHWs in the North Pacific. Given the importance of the horizontal transport of DIC in mitigating the impact of PMHWs on the uptake of CO 2 , we propose that studies should address how PMHWs and ocean circulation are interconnected in this region.

Discussion
We show that PMHWs (> 38 days and > 2.3°C anomalies as defined in this study) most frequently occur in oceanic regions of major importance for the global carbon cycle: the Tropical Pacific CO 2 source area, and the CO 2 sink regions of the North Pacific, the North Atlantic and the mid-high latitude Southern Ocean. However, over the study period 1985-2017, PMHWs have impacted air-sea CO 2 exchange only in the North and tropical Pacific CO 2 sensitive areas. The processes of this interplay are provided in the schematic of Fig. 3.
In the North Pacific CO 2 sink, PMHW events cause a reduction in the CO 2 sink as a result of the net effect of two competing mechanisms: extreme higher-than-average temperature and anomalous DIC advection. The former causes a reduction in the solubility of CO 2 in ocean water, thereby reducing the ocean uptake of CO 2 whereas the latter increases the ocean CO 2 uptake-through decreased levels of DIC driven by horizontal advection-and as such attenuates the impact of the thermal effect. Overall, the thermal effect dominates the advection effect in our study, leading to a net reduction in the air-to-sea CO 2 flux density during a PMHW event of about 29 ± 11% (Fig. 3a).
In the Tropical Pacific where CO 2 release takes place, the CO 2 outgassing is significantly attenuated during PMHWs with a reduction in CO 2 release from the ocean to the atmosphere of about 40 ± 9%. In this region, PMHWs are associated with ENSO 5 and previous studies have investigated the mechanisms explaining this change, which is mainly driven by a change in vertical ocean circulation [28][29][30][31][32] . During PMHWs (Fig. 3b), eastward propagating Kelvin waves that depress the thermocline in the east together with a concurrent weakening of easterly winds, and the extension of the western Pacific warm pool to the east, reduce the upwelling of DIC leading to a net decrease in the sea-to-air CO 2 flux density. We have only considered the average sea-to-air CO 2 flux density response to PMHWs in the North and Tropical Pacific. Future studies should investigate the spatial variability of the response in these two regions. The results in the North Pacific CO 2 sink complete previous studies of ENSO-related PMHWs in the Tropical Pacific, and together with the new results obtained in this study thus provide a comprehensive view on the interplay between PMHWs and CO 2 sensitive areas in the North and Tropical Pacific as illustrated in Fig. 3.
Over the 1985-2017 period, considering the duration and affected areas in the North Pacific and Tropical Pacific, we derive an integrated flux anomaly linked to PMHWs of 13 ± 6 TgC (Teragrams of Carbon = 10 12 grams of Carbon) and −118 ± 25 TgC in the North Pacific and Tropical Pacific respectively. The climatological integrated fluxes expected during the same time and within the same area are −119 ± 12 TgC and 307 ± 31 TgC respectively. This corresponds to a flux anomaly of 11 ± 5% and 39 ± 11% for the respective regions. While the effect in absolute terms appears small at first sight, particularly in the North Pacific (only about 0.5% of all annual marine uptake of anthropogenic CO 2 37 ), increasing duration and intensity of the heat waves become relevant for closing the current carbon budget imbalance (currently in the order of ±100-500 TgC/yr for the past decades and recent years (see Table 6 in Friedlingstein et al. 37 ) and should therefore receive attention in future budgets.
In this study, we have introduced the concept of intense and longlasting MHWs, here defined as PMHWs, with a criterion based on the duration and the mean SST anomaly (i.e., duration > 38 days and mean SST anomaly > 2.3 degrees Celsius). Our results have shown that PMHWs occur predominantly in the major CO 2 source and sink areas. Shorter-lasting MHW events are more widespread than PMHWs 10,19 and further ensemble-based analyses are needed to understand on how their impact might affect air-sea CO 2 exchanges at global scale. Except for one single product (JENA), all other observational-based products are available at monthly resolution which makes such investigation not possible. Future analyses based on an ensemble of CO 2 products of high temporal resolution (i.e. < 30 day) would be needed to study the integrated impact of shorter extreme events on ocean CO 2 fluxes.
In the North Atlantic and the mid-high latitude Southern Ocean CO 2 sinks, no change is observed for the CO 2 fluxes due to PMHWs, suggesting that the thermal effect and non-thermal effects cancel each other out. To further analyze this hypothesis, the use of reanalyzes is needed. Currently, long-term in-situ observations for pCO 2 , temperature, salinity, DIC and ALK are lacking in these areas preventing the establishment of robust air-sea CO 2 flux products, climatologies and reanalysis skill validation for such a process study. Hence, PMHWinduced impacts on CO 2 fluxes mechanisms in the Atlantic and midhigh latitude Southern Ocean remain unanswered, and increased monitoring efforts are needed accordingly.
There has been a statistically significant increase (P = .08) in the attenuation in F CO2 due to PMHWs between the 1985-1995 and the 2007-2017 period in the North Pacific, and we have compared the distribution of the ensemble of the 4 observation-based products. On the contrary, there has been no statistically significant change in the Tropical Pacific (Fig. 4). Similarly, PMHWs have also increased in intensity in the North Pacific over the last decades (linear trend= 0.11 ± 0.03°C/decade, P = .002), while their strength remains similar in the Tropics (linear trend= 0.04 ± 0.09°C/decade, P = .66) (Fig. S4). Based on the results of our process study, we can develop the hypothesis that the reported increase in the intensity of PMHWs has potentially amplified the outgassing of CO 2 over the 1985-2017 period, and that the competing mechanism, i.e., anomalous advection of DIC, was unable to counter-interact the thermal effect over this time scale. However, we cannot test such hypothesis as it would require a decomposition of F CO2 and DIC budgets from 1985 to 2017, the latter being currently not estimated by the reanalysis over this period. MHWs are projected to become stronger, more frequent and longer lasting in a warming climate 1,9,12,13 . Therefore, it is crucial to understand how PMHWs and F CO2 interact over longer time scales if we want to further unravel the evolution of the oceanic carbon cycle under climate change.

Observation-based products of F CO2
In this study, we use an ensemble of four observation-based products to quantify the impact of PMHWs on F CO2 in the three largest oceanic CO 2 sink and the largest oceanic CO 2 source in the Tropical Pacific from 1985 to 2017. Here, we provide a brief outline of the chosen products. More detail can be found in their respective publications.
The first observation-based product, from the Max Plank Institute for Meteorology (hereinafter denoted MPI) 47,48 , is based on a selforganizing map-feed-forward network that reconstructs the sea surface partial pressure of CO 2 (spCO 2 ) from various environmental predictor data. In a first step, the ocean is divided into biogeochemical regions of similar spCO 2 properties (making use of a spCO 2 climatology) and in a second step the non-linear relationship between auxiliary driver data and sparse observations is reconstructed to fill measurement gaps. The period of analysis is from 1982 to 2019 at monthly intervals and with a spatial resolution of 1°× 1°. It is based on a collection of ship and mooring spCO 2 measurements assembled by the Surface Ocean CO 2 Atlas (SOCAT) version 2020 [49][50][51][52] .
The second observation-based product, from Copernicus Marine Environmental Monitoring Service (hereinafter denoted CMEMS) 36 , is from an ensemble-based forward feed neural network that reconstruct change in spCO 2 from environmental predictor data. The period of analysis is from 1985 to 2018 at monthly intervals and with a spatial resolution of 1°× 1°. It is based on a collection of ship and mooring spCO 2 measurements assembled by the Surface Ocean CO 2 Atlas (SOCAT) version 2019 [49][50][51][52] .
The third observation-based product, from the Council for Scientific and Industrial Research (hereinafter denoted CSIR) 53 , is from a machine-learning ensemble average of six two-step clustering-regression models that reconstruct change in spCO 2 from environmental predictor data. The period of analysis is from 1982 to 2019 at monthly intervals and with a spatial resolution of 1°× 1°. It is based on a collection of ship and mooring spCO2 measurements assembled by the Surface Ocean CO 2 Atlas (SOCAT) version 2019 [49][50][51][52] .
The fourth observation-based product, from the Max Plank Institute for Biogeochemistry (hereinafter denoted JENA) 54 , is from an observation-driven ocean mixed-layer scheme that reconstruct change in spCO 2 by fitting a data-driven diagnostic model of ocean mixedlayer biogeochemistry to surface-ocean CO 2 partial pressure data from the SOCAT version 2019 [49][50][51][52] . The period of analysis is from 1957 to 2019 at daily intervals and with a spatial resolution of 2.5°× 2°. The daily fields were averaged into monthly fields.
Finally, to evaluate the skill of the BGC reanalysis in estimating F CO2 anomalies associated with PMHWs in the North Pacific CO 2 sink, we use an additional observation-based product, from the Japan Meteorology Agency (hereinafter denoted JMA) 55 , which is excluded from the spCO 2 ensemble as its period of analysis is shorter than the previously listed products, i.e. from 1990 to 2018. This product is based on multiple linear regressions that reconstruct change in spCO 2 from a set of environmental drivers. The temporal resolution is monthly intervals and with a spatial resolution of 1°× 1°. It is based on a collection of ship and mooring spCO 2 measurements assembled by the Surface Ocean CO 2 Atlas (SOCAT) version 2019 [49][50][51][52] .
Estimates of the air-to-sea flux density of CO 2 from spCO 2 data In the five observation-based products and the BGC reanalysis, the airto-sea CO 2 flux density (F CO2 ) is generated from spCO 2 data using the gas exchange formulation 56 , where α is the CO 2 solubility in seawater, k, a gas transfer coefficient, pCO 2atm is the atmospheric partial pressure of CO 2 and spCO 2 is the sea surface partial pressure of CO 2 . Here, negative values of F CO2 indicate uptake of CO 2 from the atmosphere to the ocean, while positive values indicate outgassing of CO 2 from the ocean to the atmosphere. Each product performs their own calculation of the fluxes and the methods are described in the respective publications.

Calculation of 2009-2017 monthly anomalies
In the reanalysis and the observation-based products, monthly anomalies (hereinafter denoted with a prime) are computed by removing a climatological value (hereinafter denoted with an overbar). The climatological value corresponds to the sum of a long-term linear trend and a monthly mean value. The monthly mean values are computed from the detrended monthly data. Calculation of 1985-2017 percent F CO2 anomalies The percent F CO2 anomalies during PMHWs and for the 1985-2017 period correspond to the monthly F CO2 anomalies divided by the monthly F CO2 climatological values. The anomalies and climatological values were computed following the method detailed previously, with the exception that monthly mean values were only computed from the detrended monthly data over the 1985-1995 period. During this decade, the number of PMHWs per year, and near globally, were the lowest of the 1985-2017 period (see Fig. S5). By calculating the anomalies relative to this "reference" decade, we make sure that the percent F CO2 anomalies represent a change with respect to oceanic conditions not impacted by PMHWs. In Fig. 1b, we represent, in each CO 2 sink/source, an ensemble of four 1985-2017 trimmed mean percent F CO2 anomalies derived from the observation-based products. We use the trimmed mean instead of the mean because it is a robust estimator of central tendency and provides a better estimation of the location of the bulk of the data than the mean when the distribution is asymmetric, which is the case here. More precisely, we use a 5% trimmed mean, i.e., the lowest 5% and the highest 5% of the data are excluded.

Taylor expansion of F CO2 anomalies
To determine the driving mechanisms causing F CO2 anomalies during PMHWs in the North Pacific CO 2 sink, we calculate a first-order Taylor series expansion of F CO2 anomalies in terms of its driving parameters (i.e., wind, upper ocean temperature, salinity, dissolved inorganic carbon (DIC), alkalinity (ALK) and the atmospheric partial pressure of CO 2 ) 28,44 .
First, we performed the linear Taylor decomposition of Eq. (1): The right-hand-side terms represent the contribution to F CO2 ' of gas transfer and solubility anomalies, atmospheric pCO 2 anomalies and spCO 2 anomalies. Note that the temperature dependence of k and α cancel each other, and (kα)′ is mainly driven by variations in wind speed 28,44 .
The spCO 2 anomalies are further decomposed into contributions from sea surface temperature anomalies (SST'), sea surface dissolved inorganic carbon anomalies (SDIC'), sea surface alkalinity anomalies (SALK') and sea surface salinity anomalies (SSS'), neglecting the second-order terms 28,44,45,57 : Substituting Eq. (3) into Eq. (2) gives the contributions of all parameters to F CO2 ' in a single expression: Here, the sea surface quantities correspond to the quantities at the first level of the ocean reanalysis estimates (z~−0.50 m). Following Doney et al. 28 , the partial derivatives in Eq. (4) were computed off-line at each grid point, taking SDIC as an example, as: where spCO2 values are calculated using the seacarb program for R (https://CRAN.R-project.org/package=seacarb).

DIC anomalies budget
In our study, we show that DIC anomalies play a significant role in controlling F CO2 anomalies during PMHWs. We therefore conduct a DIC anomalies budget to elucidate what processes controlled DIC anomalies.
In the ocean reanalysis, the changes in DIC concentration with time are described by the following equation: where ADV H and ADV _z are the horizontal and vertical advection of DIC respectively, DIFF z is the vertical diffusion of DIC, SBC are the freshwater fluxes that dilute or concentrate DIC, F CO2 is the air-sea CO 2 flux density, B is the biological activity that consumes or releases DIC (see details in Aumont et al. 58 ), and r is the climatological damping (see supplementary information). Positive values result in a net increase in DIC. All terms were computed online on a daily basis and stored for monthly averages. The DIC tendency (rate of change or trend) equation (Eq. 6) is expressed as a function of monthly anomalies and averaged over the average mixing layer observed during PMHWs in the reanalysis (indicated by angle brackets), i.e. from the surface to h4 7 m: Satellite sea surface temperature and marine heatwaves detection MHWs locations, dates of onset and durations were derived from the global daily remotely sensed National Ocean Atmospheric Administration (NOAA) Optimum Interpolation sea surface temperature V2, ¼°gridded data over 1982-2017 33,34 . This dataset is derived from the advanced very high-resolution radiometer (AVHRR). We apply a standard MHW detection algorithm 2 to the gridded SST data. More specifically, a warm event is considered as a MHW if it lasts for 5 or more days, with sea surface temperatures warmer than the 90th percentile based on a 1983-2012 historical climatology. The MHW detection algorithm is usually not performed on grid cells with periods of ice coverage longer than 5 days 10 . We therefore restrict our analysis to the area between 60°S and 60°N. For each MHW detected, the date of onset, duration and mean sea surface temperature anomaly are estimated by the MHW detection algorithm.

Calculation of anomalies associated with PMHWs
For each PMWH detected, the monthly anomalies were extracted at the model or observation-based products grid-point the closest to the PMHW location and for the entire duration of the PMHW. Then, to match the temporal resolution of the PMHW, the extracted anomalies were resampled from monthly to daily frequency through linear interpolation. The interpolated values were then averaged over the duration of the PMHW to give a single value, consistently with the other metrics derived from the MHW detection algorithm.